{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "name": "YeastNetwork.ipynb",
      "version": "0.3.2",
      "provenance": [
        {
          "file_id": "1TEidhGyqU7bsBLnRqmmwZ5161csK-k8g",
          "timestamp": 1547245813681
        }
      ],
      "collapsed_sections": [],
      "last_runtime": {
        "build_target": "",
        "kind": "local"
      }
    }
  },
  "cells": [
    {
      "metadata": {
        "id": "DhWmZgAVwDGu",
        "colab_type": "text"
      },
      "cell_type": "markdown",
      "source": [
        "Copyright 2019 Google LLC\n",
        "\n",
        "Licensed under the Apache License, Version 2.0 (the \"License\");\n",
        "you may not use this file except in compliance with the License.\n",
        "You may obtain a copy of the License at\n",
        "\n",
        "    https://www.apache.org/licenses/LICENSE-2.0\n",
        "\n",
        "Unless required by applicable law or agreed to in writing, software\n",
        "distributed under the License is distributed on an \"AS IS\" BASIS,\n",
        "WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.\n",
        "See the License for the specific language governing permissions and\n",
        "limitations under the License."
      ]
    },
    {
      "metadata": {
        "id": "Z7ERCe5N-RXk",
        "colab_type": "code",
        "colab": {}
      },
      "cell_type": "code",
      "source": [
        "import os\n",
        "import pandas as pd\n",
        "import numpy as np\n",
        "import scipy.stats\n",
        "import sklearn.linear_model"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "metadata": {
        "id": "qNt-xbVx-yuK",
        "colab_type": "code",
        "colab": {}
      },
      "cell_type": "code",
      "source": [
        "# Replace this with the location of the downloaded data and models\n",
        "datadir = '/path/to/datafiles'"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "metadata": {
        "id": "Y9WvQ7hYp99d",
        "colab_type": "text"
      },
      "cell_type": "markdown",
      "source": [
        "# **Common functions**"
      ]
    },
    {
      "metadata": {
        "id": "B5KdexHTfkSh",
        "colab_type": "text"
      },
      "cell_type": "markdown",
      "source": [
        "## **Auxiliary function for design matrices.**"
      ]
    },
    {
      "metadata": {
        "id": "QD1y3SXAI-3Q",
        "colab_type": "code",
        "colab": {}
      },
      "cell_type": "code",
      "source": [
        "def construct_difference_operators(times, only_x0=False):\n",
        "  \"\"\"Constuct matrix operators for differencing.\n",
        "\n",
        "  The times arg is the concatention of strictly increasing sequences, each\n",
        "  starting from zero.  The output operators are all len(times) X len(times)\n",
        "  matrices that operate on vectors of gene expression values that are len(times)\n",
        "  long.\n",
        "\n",
        "  Args:\n",
        "    times: np.array of time values\n",
        "    only_x0: if true, return the x0_operator, None, None.\n",
        "\n",
        "  Returns:\n",
        "    x0: operator that returns the time zero point for the timecourse\n",
        "    diff: operator that returns the finite-difference for the timescourse\n",
        "    diff_inv: operator that returns the timecourse from the finite differences\n",
        "\n",
        "  Raises:\n",
        "    ValueError:  First time point isn't zero or timeseries is garbled.\n",
        "  \"\"\"\n",
        "  num_times = len(times)\n",
        "  # Tag the zeros\n",
        "  zero_f = times == 0.0\n",
        "  if not zero_f[0]:\n",
        "    # First time point is not zero, something is wrong.\n",
        "    raise ValueError('First time point is non-zero, exiting.')\n",
        "  else:\n",
        "    old_t = -1\n",
        "    for t in times:\n",
        "      if t != 0 and t <= old_t:\n",
        "        raise ValueError('Timeseries not sorted properly.')\n",
        "      old_t = t\n",
        "\n",
        "  # The x0 operator acting on an expression vector returns a vector of the same\n",
        "  # length containing the time zero point replicated for the full timecourse,\n",
        "  # for each timecourse concatenated.\n",
        "  x0 = np.zeros((num_times, num_times))\n",
        "  for i, z in enumerate(zero_f):\n",
        "    if z:\n",
        "      col = i\n",
        "    x0[i, col] = 1.0\n",
        "  if only_x0:\n",
        "    return x0, None, None\n",
        "\n",
        "  # Shift times backward, add a zero at the end\n",
        "  shifted_times = np.append(times[1:], 0.0)\n",
        "  # Time deltas\n",
        "  ddt = 0.5 / (shifted_times - times)\n",
        "  # Tag the endpoints\n",
        "  last_f = ddt < 0.0\n",
        "  # The negative values correspond to the \"boundaries\" between experiments.\n",
        "  # Set these to zero.\n",
        "  ddt[last_f] = 0.0\n",
        "  # Prepend with a zero so that the first experiment can see \"backward\" to\n",
        "  # a previous dummy experiment.\n",
        "  ddt = np.append(0.0, ddt)\n",
        "  # This is the difference operator\n",
        "  diff = np.zeros((num_times, num_times))\n",
        "  for i in range(num_times):\n",
        "    # Typical case along the diagonal.\n",
        "    # If the ddt values are equal, this is 0.\n",
        "    # This corresponds to the usual symmetric difference.\n",
        "    diff[i, i] = ddt[i] - ddt[i+1]\n",
        "    # Backward derivative.\n",
        "    if i - 1 >= 0:\n",
        "      diff[i, i-1] = - ddt[i]\n",
        "    # Forward derivative.\n",
        "    if i + 1 < num_times:\n",
        "      diff[i, i+1] = ddt[i+1]\n",
        "  # This is the inverse of the difference operator\n",
        "  diff_inv = np.zeros((num_times, num_times))\n",
        "  # It is not invertible.  However, removing the first column\n",
        "  # and last row from each timecourse makes it invertible.\n",
        "  # We do that because we don't care about the zero point of the sequence,\n",
        "  # thus filtering columns for non-zero time points, and we don't care about\n",
        "  # the final derivative, thus filtering rows for \"not last\".\n",
        "  d_inv = np.linalg.inv(diff[np.ix_(~last_f, ~zero_f)])\n",
        "  # Fill in the values.\n",
        "  # Note that in constructing the data from the derivatives,\n",
        "  # we do not try to reconstruct x0, thus the row filter on non-zero times\n",
        "  # and we do not use the final derivative of the timecourse due to linear\n",
        "  # dependence, thus the column filter on \"not last\".\n",
        "  # diff_inv . diff . X == X - X0\n",
        "  diff_inv[np.ix_(~zero_f, ~last_f)] = d_inv\n",
        "\n",
        "  # Return the three operators.\n",
        "  return x0, diff, diff_inv"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "metadata": {
        "id": "s0TqinOcJ3O4",
        "colab_type": "code",
        "colab": {}
      },
      "cell_type": "code",
      "source": [
        "def add_gene_block(df_base, df_add, block_name, genes):\n",
        "  # Column-wise bind df_base (dataframe) and df_add (may be np.array).\n",
        "  # Make a multindex of (blockname, gene) where gene is taken from genes.\n",
        "  df_add.columns = pd.MultiIndex.from_tuples(\n",
        "      [(block_name, x) for x in genes])\n",
        "  return pd.concat([df_base, df_add], axis=1)"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "metadata": {
        "id": "E6OqKCzFHtfA",
        "colab_type": "code",
        "colab": {}
      },
      "cell_type": "code",
      "source": [
        "TIMECOURSE_IDX = ['TF', 'strain', 'date', 'restriction', 'mechanism']\n",
        "CHIP_IDX = TIMECOURSE_IDX + ['time']\n",
        "TIMESERIES_IDX = ['GeneName', 'SystematicName'] + CHIP_IDX\n",
        "DESIGN_IDX = CHIP_IDX + ['seq', 'rseq']"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "metadata": {
        "id": "4kVCZJ_Vh3Dg",
        "colab_type": "text"
      },
      "cell_type": "markdown",
      "source": [
        "## **Auxiliary functions for noise model**"
      ]
    },
    {
      "metadata": {
        "id": "R-1BYVLsIN5g",
        "colab_type": "code",
        "colab": {}
      },
      "cell_type": "code",
      "source": [
        "def normalize_to_t0(df_dm, x0_op):\n",
        "  \"\"\"Constructs variables relative to t0.\n",
        "\n",
        "  Args:\n",
        "    df_dm: dataframe of the design_matrix\n",
        "    x0_op: operator to produce the t=0 point for the timeseries\n",
        "\n",
        "  Returns:\n",
        "    design matrix dataframe with ratio block.\n",
        "  \"\"\"\n",
        "  genes = df_dm['r_g_median'].columns\n",
        "\n",
        "  # The basic stat used is the median.\n",
        "  # Here we compare with the zero time\n",
        "  df_dm = add_gene_block(\n",
        "      df_dm,\n",
        "      np.log(df_dm['r_g_median'] / np.matmul(x0_op, df_dm['r_g_median'])),\n",
        "      'log_ratio',\n",
        "      genes)\n",
        "\n",
        "  return df_dm"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "metadata": {
        "id": "PvGr0nRVIawy",
        "colab_type": "code",
        "colab": {}
      },
      "cell_type": "code",
      "source": [
        "def repair_crosstalk_outliers(df, x0_op):\n",
        "  \"\"\"Corrects for timeseries with severe red - green crosstalk.\n",
        "\n",
        "  Args:\n",
        "    df: dataframe of the design_matrix\n",
        "    x0_op: operator to produce the t=0 point for the timeseries\n",
        "\n",
        "  Returns:\n",
        "    design matrix dataframe with corrected ratio.\n",
        "  \"\"\"\n",
        "  # Pull out the zero points of the green channel, convert back to DataFrame.\n",
        "  green0 = pd.DataFrame(\n",
        "      np.matmul(x0_op, df['green_median']),\n",
        "      index=df.index,\n",
        "      columns=df['green_median'].columns)\n",
        "\n",
        "  # Removing the t0 point, these are the 30% quantiles of each\n",
        "  # timecourse + gene pair.  This covers 2/7 of the remaining points.\n",
        "  # The stack() turns the gene columns into rows.\n",
        "  gratio = ((df['green_median'] / green0).query('time > 0').\n",
        "            groupby(level=['TF', 'strain']).quantile(0.3).stack())\n",
        "  lratio = (df['log_ratio'].query('time > 0').\n",
        "            groupby(level=['TF', 'strain']).quantile(0.3).stack())\n",
        "\n",
        "  # Flag the genes that both have green going way up (8x) and also the\n",
        "  # uncorrected ratio going up 2x beyond that as \"crosstalk\" genes.\n",
        "  # When iterated on, this index is (TF, strain, gene) tuples.\n",
        "  crosstalk_genes = gratio[(gratio > 8.0) & (lratio > np.log(2.0))].index\n",
        "  for tf, strain, gene in crosstalk_genes:\n",
        "    print(\n",
        "        'Repairing %s, %s timecourse for gene %s.' % (tf, strain, gene))\n",
        "    tc_f = ((df.index.get_level_values('TF') == tf) &\n",
        "            (df.index.get_level_values('strain') == strain))\n",
        "    # This adds log(green / green0), which yields\n",
        "    # log_ratio = log(red) - log(green0)\n",
        "    # (the green0 is propagated across the timecourse).\n",
        "    df.loc[tc_f, ('log_ratio', gene)] += np.log(\n",
        "        df.loc[tc_f, ('green_median', gene)] / green0.loc[tc_f, gene])\n",
        "\n",
        "  return df"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "metadata": {
        "id": "T9nQ26BmKglA",
        "colab_type": "code",
        "colab": {}
      },
      "cell_type": "code",
      "source": [
        "# The typical timecourse was 0, 5, 10, 15, 20, 30, 45, 90\n",
        "TIME_POINTS = (0, 6, 12, 17.5, 25, 35, 85, 1000)\n",
        "TIME_BINS = zip(TIME_POINTS[:-1], TIME_POINTS[1:])"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "metadata": {
        "id": "uSXH-R5WIn7n",
        "colab_type": "code",
        "colab": {}
      },
      "cell_type": "code",
      "source": [
        "def clean_time_trends(df_dm):\n",
        "  \"\"\"Corrects for median timeseries behavior.\n",
        "\n",
        "  Args:\n",
        "    df_dm: dataframe of the design_matrix\n",
        "\n",
        "  Returns:\n",
        "    design matrix dataframe with new blocks.\n",
        "  \"\"\"\n",
        "  genes = df_dm['log_ratio'].columns\n",
        "\n",
        "  dfr = df_dm['log_ratio']\n",
        "  dfrc = dfr.copy()\n",
        "  times = dfr.index.get_level_values('time')\n",
        "  gev = dfr.index.get_level_values('mechanism') == 'GEV'\n",
        "  for mn, mx in TIME_BINS:\n",
        "    for g in (gev, ~gev):\n",
        "      filt = (times > mn) & (times <= mx) & g\n",
        "      dfrc[filt] -= dfr[filt].median()\n",
        "\n",
        "  df_dm = add_gene_block(\n",
        "      df_dm,\n",
        "      dfrc,\n",
        "      'log_cleaned_ratio',\n",
        "      genes)\n",
        "\n",
        "  return df_dm"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "metadata": {
        "id": "EhST0tf-rPVk",
        "colab_type": "code",
        "colab": {}
      },
      "cell_type": "code",
      "source": [
        "def build_cleaned_table(df_dm):\n",
        "  \"\"\"Augment design table and write to disk.\n",
        "\n",
        "  Args:\n",
        "    df_dm: design table dataframe\n",
        "\n",
        "  Returns:\n",
        "    design table dataframe with normalized, cleaned, and differentiated columns\n",
        "  \"\"\"\n",
        "  x0_op, _, _ = construct_difference_operators(\n",
        "      np.array(df_dm.index.get_level_values('time')), only_x0=True)\n",
        "\n",
        "  df_dm = normalize_to_t0(df_dm, x0_op)\n",
        "\n",
        "  df_dm = repair_crosstalk_outliers(df_dm, x0_op)\n",
        "\n",
        "  df_dm = clean_time_trends(df_dm)\n",
        "\n",
        "  df_dm = construct_noise_model(df_dm)\n",
        "\n",
        "  return df_dm"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "metadata": {
        "id": "G7-4Xv9TK_yA",
        "colab_type": "code",
        "colab": {}
      },
      "cell_type": "code",
      "source": [
        "# For the normal distribution, the central quintile is +/- 0.253, approximately.\n",
        "# We scale to infer the standard deviation.\n",
        "CENTRAL_QUINTILE_SCALE = 1.0 / (scipy.stats.norm.ppf(0.6) -\n",
        "                                scipy.stats.norm.ppf(0.4))\n",
        "\n",
        "def central_quintile_error(dfc, axis=0):\n",
        "  # Use central quintile to define the noise level.\n",
        "  return (dfc.quantile(0.6, axis=axis) -\n",
        "          dfc.quantile(0.4, axis=axis)) * CENTRAL_QUINTILE_SCALE\n",
        "\n",
        "\n",
        "def log_cleaned_gene_error(df):\n",
        "  return central_quintile_error(df.query('time > 0')['log_cleaned_ratio'])\n",
        "\n",
        "\n",
        "def chip_residual_error(df):\n",
        "  return central_quintile_error(df[df != 0], axis=1)"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "metadata": {
        "id": "FR2qlqAeNF0n",
        "colab_type": "code",
        "colab": {}
      },
      "cell_type": "code",
      "source": [
        "def apply_soft_thresholding(df, noise):\n",
        "  return np.median(\n",
        "      np.array([np.array(df + noise),\n",
        "                np.array(df - noise),\n",
        "                np.zeros(df.shape)]),\n",
        "      axis=0)\n",
        "\n",
        "\n",
        "def apply_hard_thresholding(df, noise):\n",
        "  hth = df.copy()\n",
        "  hth[np.abs(df) < noise] = 0.0\n",
        "  return hth"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "metadata": {
        "id": "BfGLwnMuN10G",
        "colab_type": "code",
        "colab": {}
      },
      "cell_type": "code",
      "source": [
        "def count_timecourse_detections(df):\n",
        "  # Select active genes for each timecourse.\n",
        "  # A gene is active in a timecourse if it ever passes hard thresholding.\n",
        "  # Drop the unused names from the initial index.\n",
        "  drop_names = [x for x in df.index.names if x not in TIMECOURSE_IDX]\n",
        "\n",
        "  return (\n",
        "      # Significant data points (bool)\n",
        "      ((df != 0).\n",
        "       # Expanded index\n",
        "       reset_index().\n",
        "       # Group by timecourses\n",
        "       groupby(TIMECOURSE_IDX).\n",
        "       # Sum just counts significant points\n",
        "       sum()).\n",
        "      # Drop unwanted terms in expanded index\n",
        "      drop(drop_names, axis=1, errors='ignore'))\n",
        "\n",
        "\n",
        "def flag_non_zero_timecourses(df):\n",
        "  return count_timecourse_detections(df) != 0\n",
        "\n",
        "\n",
        "def count_consecutive_timecourse_detections(df):\n",
        "  # Check for consecutive points passing thresholding.\n",
        "  # Shift forward\n",
        "  shifted = np.zeros_like(df)\n",
        "  shifted[:-1] = (df[1:] != 0)\n",
        "  # consec == 0 for first and last points, automatically.\n",
        "  consec = (df != 0) * shifted\n",
        "  return count_timecourse_detections(consec)\n",
        "\n",
        "\n",
        "def strip_timecourses(df, sig_tc, significant=True):\n",
        "  \"\"\"Strip timecourses according to a set of significant ones.\n",
        "\n",
        "  Args:\n",
        "    df: design matrix dataframe to be stripped.\n",
        "    sig_tc: boolean design with timecourse index and gene columns.\n",
        "    significant: if true, strip the significant timecourses,\n",
        "                 otherwise, strip the insignificant ones.\n",
        "\n",
        "  Returns:\n",
        "    dataframe of stripped design matrix\n",
        "  \"\"\"\n",
        "  dfn = df.copy()\n",
        "  for idx, sig in sig_tc.iterrows():\n",
        "    # This builds the query of index_name==index_value, joined with ' and '.\n",
        "    query = ' and '.join(['%s==\"%s\"' % t for t in zip(TIMECOURSE_IDX, idx)])\n",
        "    filt = dfn.query(query)\n",
        "    dfn.loc[filt.index] *= np.array(~sig if significant else sig, dtype=float)\n",
        "  return dfn"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "metadata": {
        "id": "Q8QQ8_D3Jnhv",
        "colab_type": "code",
        "colab": {}
      },
      "cell_type": "code",
      "source": [
        "MINIMUM_THRESHOLDING = 0.1\n",
        "CHIP_WEIGHT = 0.2\n",
        "def construct_noise_model(df_dm):\n",
        "  \"\"\"Builds a gene X chip noise model.\n",
        "\n",
        "  Args:\n",
        "    df_dm: dataframe of the design_matrix\n",
        "\n",
        "  Returns:\n",
        "    design matrix dataframe with new blocks.\n",
        "  \"\"\"\n",
        "  genes = df_dm['log_cleaned_ratio'].columns\n",
        "  gene_error = log_cleaned_gene_error(df_dm)\n",
        "\n",
        "  # Thresholding factor.\n",
        "  threshold_factor = np.sqrt(2.0 * np.log(\n",
        "      len(df_dm.query('time > 0'))))\n",
        "  threshold_error = np.maximum(\n",
        "      gene_error * threshold_factor,\n",
        "      MINIMUM_THRESHOLDING)\n",
        "\n",
        "  hth = apply_hard_thresholding(\n",
        "      df_dm['log_cleaned_ratio'], threshold_error)\n",
        "  # For each timecourse, zero out the significant genes.\n",
        "  sig_tc = flag_non_zero_timecourses(hth)\n",
        "  dfn = strip_timecourses(\n",
        "      df_dm['log_cleaned_ratio'], sig_tc)\n",
        "  # Compute chip-level error via central quintile.\n",
        "  chip_error = chip_residual_error(dfn)\n",
        "  filt = chip_error.index.get_level_values('time') == 0.0\n",
        "  # Fill in 0.0 for the t=0 points. t=0 will thus have gene-level noise\n",
        "  chip_error[filt] = 0.0\n",
        "\n",
        "  noise_df = pd.DataFrame(\n",
        "      np.sqrt(np.ones(dfn.shape) * gene_error.values.reshape(1, -1)**2 +\n",
        "              CHIP_WEIGHT *\n",
        "              np.ones(dfn.shape) * chip_error.values.reshape(-1, 1)**2),\n",
        "      index=chip_error.index,\n",
        "      columns=gene_error.index)\n",
        "\n",
        "  # Normalize the non-zero points properly.\n",
        "  nonzero = noise_df.query('seq != 0')\n",
        "  noise_df.loc[nonzero.index] /= np.sqrt(1.0 + CHIP_WEIGHT)\n",
        "\n",
        "  # Output the noise model.\n",
        "  df_dm = add_gene_block(\n",
        "      df_dm,\n",
        "      noise_df,\n",
        "      'log_noise_model',\n",
        "      genes)\n",
        "\n",
        "  return df_dm"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "metadata": {
        "id": "MnbdIrmvMQ_c",
        "colab_type": "code",
        "colab": {}
      },
      "cell_type": "code",
      "source": [
        "LOG_SUFFIXES = ['noise_model',\n",
        "                'ratio',\n",
        "                'cleaned_ratio',\n",
        "                'cleaned_ratio_zth',\n",
        "                'cleaned_ratio_hth',\n",
        "                'cleaned_ratio_th',\n",
        "                'cleaned_ratio_zth2d',\n",
        "                'cleaned_ratio_hth2d',\n",
        "                'cleaned_ratio_th2d',\n",
        "                'cleaned_ratio_zth2dfilt',\n",
        "                'cleaned_ratio_hth2dfilt',\n",
        "                'cleaned_ratio_th2dfilt']"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "metadata": {
        "id": "73UJYJ1DqiW0",
        "colab_type": "text"
      },
      "cell_type": "markdown",
      "source": [
        "# **Process \"raw\" data**"
      ]
    },
    {
      "metadata": {
        "id": "my9L8TeFeemv",
        "colab_type": "text"
      },
      "cell_type": "markdown",
      "source": [
        "## **Read in the raw data and verify that the red to green ratio column is correctly computed.**"
      ]
    },
    {
      "metadata": {
        "id": "muiWwESMCxIw",
        "colab_type": "code",
        "colab": {}
      },
      "cell_type": "code",
      "source": [
        "# Read in the raw data.\n",
        "\n",
        "with open(os.path.join(datadir, 'yeast_raw_data_table_20180826.tsv')) as f:\n",
        "  df_raw = pd.read_csv(f, sep='\\t')"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "metadata": {
        "id": "boDRadhBGZP1",
        "colab_type": "code",
        "colab": {}
      },
      "cell_type": "code",
      "source": [
        "# The counts are red(mutant) and green(control).  r_g_ratio is the ratio of these, clipped at values of 2.0.\n",
        "\n",
        "d = np.maximum(df_raw['rProcessedSignal'], 2.0) / np.maximum(df_raw['gProcessedSignal'], 2.0) - df_raw['r_g_ratio']\n",
        "min(d), max(d)"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "metadata": {
        "id": "1IlFonJPerar",
        "colab_type": "text"
      },
      "cell_type": "markdown",
      "source": [
        "## **Aggregate across spots on the microarrays.**"
      ]
    },
    {
      "metadata": {
        "id": "fe4R8FVGXoTX",
        "colab_type": "code",
        "colab": {}
      },
      "cell_type": "code",
      "source": [
        "# Aggregate across the multiple spots (usually 2) on the microarrays.\n",
        "\n",
        "df_ts = df_raw.pivot_table(index=TIMESERIES_IDX,\n",
        "                           values=['r_g_ratio',\n",
        "                                   'rProcessedSignal',\n",
        "                                   'gProcessedSignal'],\n",
        "                           aggfunc=[len, np.median, np.min, np.max, np.std])"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "metadata": {
        "id": "r1-F6nqKAi0_",
        "colab_type": "code",
        "colab": {}
      },
      "cell_type": "code",
      "source": [
        "# Save memory.\n",
        "del df_raw"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "metadata": {
        "id": "t0cpPjVIe_Oc",
        "colab_type": "text"
      },
      "cell_type": "markdown",
      "source": [
        "## **Trim the most egregious outliers.**"
      ]
    },
    {
      "metadata": {
        "id": "8WivoCgoMtJh",
        "colab_type": "code",
        "colab": {}
      },
      "cell_type": "code",
      "source": [
        "MAX_SPOT_RATIO = 4.0\n",
        "JUMP_THRESHOLD = 4.0\n",
        "def trim_timeseries_outliers(df_ts):\n",
        "  \"\"\"Replace outliers in timeseries with bracketing averages.\n",
        "\n",
        "  Args:\n",
        "    df_ts: dataframe with gene X experiment X time per row\n",
        "\n",
        "  Returns:\n",
        "    dataframe where timeseries outliers have been replaced by\n",
        "    bracketing averages.\n",
        "  \"\"\"\n",
        "  # Shift times down by one to get \"next\" time.\n",
        "  next_t = pd.Series(\n",
        "      np.append(np.array(df_ts.iloc[1:].index.get_level_values('time')), 0),\n",
        "      index=df_ts.index)\n",
        "  # Values shifted forward and backward by one in time.\n",
        "  next_val = pd.Series(\n",
        "      np.append(np.array(df_ts.iloc[1:][('median', 'r_g_ratio')]), 0),\n",
        "      index=df_ts.index)\n",
        "  prev_val = pd.Series(\n",
        "      np.append(0, np.array(df_ts.iloc[:-1][('median', 'r_g_ratio')])),\n",
        "      index=df_ts.index)\n",
        "  orig_val = df_ts[('median', 'r_g_ratio')]\n",
        "  geometric_mean = np.sqrt(prev_val * next_val)\n",
        "\n",
        "  # filters for first and last in a timecourse.\n",
        "  first_t = df_ts.index.get_level_values('time') == 0\n",
        "  last_t = df_ts.index.get_level_values('time') > next_t\n",
        "\n",
        "  # Filter for genes on a chip where the spots disagree.\n",
        "  filt = (df_ts[('amax', 'r_g_ratio')] /\n",
        "          np.maximum(df_ts[('amin', 'r_g_ratio')], 0.001) >\n",
        "          MAX_SPOT_RATIO)\n",
        "\n",
        "  # Replace the cross-spot medians for the inconsistent genes.\n",
        "  m_filt = filt & ~first_t & ~last_t\n",
        "  df_ts.loc[m_filt, ('median', 'r_g_ratio')] = geometric_mean[m_filt]\n",
        "\n",
        "  f_filt = filt & first_t\n",
        "  df_ts.loc[f_filt, ('median', 'r_g_ratio')] = next_val[f_filt]\n",
        "\n",
        "  l_filt = filt & last_t\n",
        "  df_ts.loc[l_filt, ('median', 'r_g_ratio')] = prev_val[l_filt]\n",
        "\n",
        "  print('%d / %d spot-discrepant outliers trimmed. ' % (\n",
        "      sum(filt), len(filt)))\n",
        "\n",
        "  # Check for outliers from the time-trend in r_g_ratio\n",
        "  m_filt = ~first_t & ~last_t\n",
        "  ratio1 = orig_val / prev_val\n",
        "  ratio2 = next_val / orig_val\n",
        "  # outlier is a non-endpoint location where there was movement\n",
        "  # up then down or down then up, each over threshold.\n",
        "  outlier_filt = m_filt & (\n",
        "      ((ratio1 > JUMP_THRESHOLD) & (ratio2 < 1. / JUMP_THRESHOLD) |\n",
        "       (ratio1 < 1. / JUMP_THRESHOLD) & (ratio2 > JUMP_THRESHOLD)))\n",
        "  print('%d / %d time-series trend outliers trimmed.' % (\n",
        "               np.sum(outlier_filt), len(outlier_filt)))\n",
        "  # Replace them with the geometric mean of their neighbors.\n",
        "  df_ts.loc[outlier_filt, ('median', 'r_g_ratio')] = geometric_mean[\n",
        "      outlier_filt]\n",
        "\n",
        "  return df_ts\n"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "metadata": {
        "id": "QclLNDNkhNXp",
        "colab_type": "code",
        "colab": {}
      },
      "cell_type": "code",
      "source": [
        "df_ts = trim_timeseries_outliers(df_ts)"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "metadata": {
        "id": "7-RrdsFGfJJJ",
        "colab_type": "text"
      },
      "cell_type": "markdown",
      "source": [
        "## **Build a design table with rows as experiment time points and columns as transcripts.**"
      ]
    },
    {
      "metadata": {
        "id": "6uxqekKrhRRu",
        "colab_type": "code",
        "colab": {}
      },
      "cell_type": "code",
      "source": [
        "def design_table(df_tsm):\n",
        "  \"\"\"Generates the design matrix for linear modeling.\n",
        "\n",
        "  Args:\n",
        "    df_tsm: dataframe of the full dataset\n",
        "\n",
        "  Returns:\n",
        "    dataframe containing design matrix\n",
        "  \"\"\"\n",
        "  df_tsmf = df_tsm[[('median', 'r_g_ratio'),\n",
        "                    ('median', 'rProcessedSignal'),\n",
        "                    ('median', 'gProcessedSignal')]]\n",
        "  df_tsmf.columns = ['r_g_median', 'red_median', 'green_median']\n",
        "\n",
        "  df_dm = pd.pivot_table(\n",
        "      df_tsmf.reset_index(),\n",
        "      index=CHIP_IDX,\n",
        "      columns=['GeneName'])\n",
        "\n",
        "  return df_dm"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "metadata": {
        "id": "emqRcY_GIGE5",
        "colab_type": "code",
        "colab": {}
      },
      "cell_type": "code",
      "source": [
        "def timeseries_sequence(times):\n",
        "  \"\"\"Compute the timecourse sequence number of the time points, and reverse.\"\"\"\n",
        "  diffs = np.append(times[1:], 0) - times\n",
        "  i = 0\n",
        "  seq = []\n",
        "  for x in diffs:\n",
        "    seq.append(i)\n",
        "    if x <= 0:\n",
        "      i = 0\n",
        "    else:\n",
        "      i += 1\n",
        "\n",
        "  # Now in reverse.\n",
        "  rtimes = times[::-1]\n",
        "  rdiffs = np.append(0, rtimes[:-1]) - rtimes\n",
        "  i = 0\n",
        "  rseq = []\n",
        "  for x in rdiffs:\n",
        "    if x <= 0:\n",
        "      i = 0\n",
        "    else:\n",
        "      i += 1\n",
        "    rseq.append(i)\n",
        "\n",
        "  # Don't forget to re-reverse rseq\n",
        "  return seq, rseq[::-1]\n"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "metadata": {
        "id": "_CkdfVWCjVkS",
        "colab_type": "code",
        "colab": {}
      },
      "cell_type": "code",
      "source": [
        "df_dm = design_table(df_ts)\n",
        "\n",
        "seq, rseq = timeseries_sequence(\n",
        "      np.array(df_dm.index.get_level_values('time')))\n",
        "\n",
        "df_dm['seq'] = seq\n",
        "df_dm['rseq'] = rseq\n",
        "\n",
        "df_dm = df_dm.reset_index().set_index(DESIGN_IDX).sort_index()"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "metadata": {
        "id": "H34wSGbHmrbH",
        "colab_type": "code",
        "colab": {}
      },
      "cell_type": "code",
      "source": [
        "df_dm.shape"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "metadata": {
        "id": "nDhpC7Kkmzhp",
        "colab_type": "code",
        "colab": {}
      },
      "cell_type": "code",
      "source": [
        "# There a few mismatches between the transcripts on GEV vs. ZEV microarrays.\n",
        "# Reduce to the common set.\n",
        "\n",
        "df_dm.dropna(axis=1, inplace=True)"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "metadata": {
        "id": "wuRgCkngrIdG",
        "colab_type": "code",
        "colab": {}
      },
      "cell_type": "code",
      "source": [
        "df_dm.shape"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "metadata": {
        "id": "WrcDmJoTB32v",
        "colab_type": "code",
        "colab": {}
      },
      "cell_type": "code",
      "source": [
        "# Save memory.\n",
        "del df_ts"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "metadata": {
        "id": "uAheqq8eiAHu",
        "colab_type": "text"
      },
      "cell_type": "markdown",
      "source": [
        "## **Build the noise model and \"cleaned\" dataset.**"
      ]
    },
    {
      "metadata": {
        "id": "7kiWdbYds-fP",
        "colab_type": "code",
        "colab": {}
      },
      "cell_type": "code",
      "source": [
        "# Data cleaning and noise model depend on the full dataset.\n",
        "# To reproduce the first round of results, we need to filter out\n",
        "# the rounds of validation experiments.\n",
        "\n",
        "validation_f = pd.to_datetime(df_dm.index.get_level_values('date')) > '2017-07-01'"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "metadata": {
        "id": "XvfBfonlKspq",
        "colab_type": "code",
        "colab": {}
      },
      "cell_type": "code",
      "source": [
        "# This will recapitulate the original training set,\n",
        "# before validation experiments were performed.\n",
        "# The outputs of cleaning and noise model are slightly different.\n",
        "#\n",
        "# df_dm = build_cleaned_table(df_dm[~validation_f])\n",
        "\n",
        "# This recapitulates the full dataset.\n",
        "df_dm = build_cleaned_table(df_dm)"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "metadata": {
        "id": "DtIcrHWnL85Q",
        "colab_type": "code",
        "colab": {}
      },
      "cell_type": "code",
      "source": [
        "df_dm.shape"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "metadata": {
        "id": "uVHFUc2miR9H",
        "colab_type": "text"
      },
      "cell_type": "markdown",
      "source": [
        "## **Apply noise model to generate \"thresholded\" dataset.**"
      ]
    },
    {
      "metadata": {
        "id": "4OB8MQygONV-",
        "colab_type": "code",
        "colab": {}
      },
      "cell_type": "code",
      "source": [
        "def apply_all_thresholding(df_dm, min_thresholding=0.1):\n",
        "  \"\"\"Adds thresholded data blocks to design matrix.\n",
        "\n",
        "  Args:\n",
        "    df_dm: dataframe of the design_matrix\n",
        "    min_thresholding: floor for amount of thresholding\n",
        "\n",
        "  Returns:\n",
        "    design matrix dataframe with new blocks.\n",
        "  \"\"\"\n",
        "  genes = df_dm['log_ratio'].columns\n",
        "  # Soft thresholding.\n",
        "  # Consider the max of N draws from a normal distribution.\n",
        "  n_eff = len(df_dm.query('time > 0'))\n",
        "  soft_threshold_factor = np.sqrt(2.0 * np.log(n_eff))\n",
        "\n",
        "  print(\n",
        "      'Thresholding at %.3f sigma for %d samples. ' % (\n",
        "      soft_threshold_factor,\n",
        "      n_eff))\n",
        "\n",
        "  gene_error = df_dm['log_noise_model'].query('time==0').median()\n",
        "\n",
        "  threshold_error = np.maximum(\n",
        "      gene_error * soft_threshold_factor,\n",
        "      min_thresholding)\n",
        "\n",
        "  df_dm = apply_hard_and_soft_thresholding(\n",
        "      df_dm, threshold_error, genes, 'th')\n",
        "\n",
        "  threshold_error_2d = np.maximum(\n",
        "      df_dm['log_noise_model'] * soft_threshold_factor,\n",
        "      min_thresholding)\n",
        "\n",
        "  df_dm = apply_hard_and_soft_thresholding(\n",
        "      df_dm, threshold_error_2d, genes, 'th2d')\n",
        "\n",
        "  df_dm = apply_hard_and_soft_thresholding(\n",
        "      df_dm, threshold_error_2d, genes, 'th2dfilt', aggressive=True)\n",
        "\n",
        "  return df_dm"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "metadata": {
        "id": "74IvBOWzPB-V",
        "colab_type": "code",
        "colab": {}
      },
      "cell_type": "code",
      "source": [
        "def apply_hard_and_soft_thresholding(df_dm, noise, genes, suffix,\n",
        "                                     aggressive=False):\n",
        "  \"\"\"Adds thresholded blocks to the design matrix.\n",
        "\n",
        "  Args:\n",
        "    df_dm: dataframe of the design_matrix\n",
        "    noise: compatible dataframe (perhaps broadcast) with noise model\n",
        "    genes: list of gene names\n",
        "    suffix: suffix for the names of the new blocks\n",
        "    aggressive: if true, apply simple timecourse-level filter as well\n",
        "\n",
        "  Returns:\n",
        "    design matrix dataframe with new blocks.\n",
        "  \"\"\"\n",
        "  # Here we hard threshold the values\n",
        "  ratio_hth = apply_hard_thresholding(df_dm['log_cleaned_ratio'], noise)\n",
        "\n",
        "  tc_counts = count_timecourse_detections(ratio_hth)\n",
        "  tc_consec_counts = count_consecutive_timecourse_detections(ratio_hth)\n",
        "  last_counts = count_timecourse_detections(ratio_hth.query('rseq==0'))\n",
        "\n",
        "  print('%d timecourses with signal (%s)' % (\n",
        "               (tc_counts > 0).sum().sum(),\n",
        "               suffix))\n",
        "  print('%d timecourses with one detection (%s)' % (\n",
        "               (tc_counts == 1).sum().sum(),\n",
        "               suffix))\n",
        "  print('%d timecourses with only final detection (%s)' % (\n",
        "               ((tc_counts == 1) & last_counts).sum().sum(),\n",
        "               suffix))\n",
        "  print('%d timecourses with no detection gaps (%s)' % (\n",
        "               (tc_consec_counts == (tc_counts - 1)).sum().sum(),\n",
        "               suffix))\n",
        "\n",
        "  if aggressive:\n",
        "    # Non-final single detection is bad.\n",
        "    notlast = (tc_counts == 1) & (last_counts == 0)\n",
        "    # Non-consecutive detection is bad.\n",
        "    nonconsec = tc_consec_counts < (tc_counts - 1)\n",
        "\n",
        "    bad_tc = notlast | nonconsec\n",
        "\n",
        "    good_tc = (tc_counts > 0) & ~bad_tc\n",
        "    print('%d timecourses with no detection gaps or singleton not last (%s)' % (\n",
        "              good_tc.sum().sum(),\n",
        "              suffix))\n",
        "\n",
        "    stripped_tc = strip_timecourses(df_dm['log_cleaned_ratio'],\n",
        "                                    good_tc,\n",
        "                                    significant=False)\n",
        "    # We need the hard thresholded data without the bad timecourses.\n",
        "    ratio_hth = apply_hard_thresholding(stripped_tc, noise)\n",
        "\n",
        "  else:\n",
        "    stripped_tc = strip_timecourses(df_dm['log_cleaned_ratio'],\n",
        "                                    tc_counts > 0,\n",
        "                                    significant=False)\n",
        "\n",
        "  df_dm = add_gene_block(\n",
        "      df_dm,\n",
        "      stripped_tc.copy(),\n",
        "      'log_cleaned_ratio_z' + suffix,\n",
        "      genes)\n",
        "\n",
        "  df_dm = add_gene_block(\n",
        "      df_dm,\n",
        "      ratio_hth,\n",
        "      'log_cleaned_ratio_h' + suffix,\n",
        "      genes)\n",
        "\n",
        "  # Here we soft threshold the values\n",
        "  ratio_th = apply_soft_thresholding(stripped_tc, noise)\n",
        "\n",
        "  df_dm = add_gene_block(\n",
        "      df_dm,\n",
        "      pd.DataFrame(ratio_th, index=df_dm.index, columns=genes),\n",
        "      'log_cleaned_ratio_' + suffix,\n",
        "      genes)\n",
        "\n",
        "  return df_dm"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "metadata": {
        "id": "TXfsCvMbt2o1",
        "colab_type": "code",
        "colab": {}
      },
      "cell_type": "code",
      "source": [
        "df_dm = apply_all_thresholding(\n",
        "        df_dm,\n",
        "        min_thresholding=MINIMUM_THRESHOLDING)"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "metadata": {
        "id": "ZizAOwptUVW6",
        "colab_type": "code",
        "colab": {}
      },
      "cell_type": "code",
      "source": [
        "df_dm['log_cleaned_ratio_th'].shape"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "metadata": {
        "id": "7Wf8OnU2ibwG",
        "colab_type": "text"
      },
      "cell_type": "markdown",
      "source": [
        "## **Generate the data table in log2 units.**"
      ]
    },
    {
      "metadata": {
        "id": "89EB08hVRMEH",
        "colab_type": "code",
        "colab": {}
      },
      "cell_type": "code",
      "source": [
        "def biologist_table(df):\n",
        "  \"\"\"Output a table in biologist format - stacked and log base 2.\"\"\"\n",
        "  base_columns = (DESIGN_IDX +\n",
        "                  ['GeneName',\n",
        "                   'green_median',\n",
        "                   'red_median',\n",
        "                   'r_g_median'])\n",
        "\n",
        "  dfs = df.stack().reset_index()\n",
        "  log_columns = ['log_' + x for x in LOG_SUFFIXES]\n",
        "  dfs[log_columns] /= np.log(2.0)\n",
        "\n",
        "  log2_columns = ['log2_' + x for x in LOG_SUFFIXES]\n",
        "  dfb = dfs[base_columns + log_columns]\n",
        "  dfb.columns = [base_columns + log2_columns]\n",
        "  return dfb"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "metadata": {
        "id": "OBmbLN3p1kLJ",
        "colab_type": "code",
        "colab": {}
      },
      "cell_type": "code",
      "source": [
        "df_bio = biologist_table(df_dm)"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "metadata": {
        "id": "WSjxNgO-1qMw",
        "colab_type": "code",
        "colab": {}
      },
      "cell_type": "code",
      "source": [
        "df_bio.shape"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "metadata": {
        "id": "A2h_yrR5qtpr",
        "colab_type": "text"
      },
      "cell_type": "markdown",
      "source": [
        "# **Read processed data**"
      ]
    },
    {
      "metadata": {
        "id": "vrGV8an3ivwm",
        "colab_type": "text"
      },
      "cell_type": "markdown",
      "source": [
        "## **Read in the processed table.**"
      ]
    },
    {
      "metadata": {
        "id": "MulgvIp11rMW",
        "colab_type": "code",
        "colab": {}
      },
      "cell_type": "code",
      "source": [
        "# Read in the final data.\n",
        "\n",
        "with open(os.path.join(datadir, 'yeast_data_table_20180826.tsv')) as f:\n",
        "  df_bio_read = pd.read_csv(f, sep='\\t')"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "metadata": {
        "id": "0kVu5KS45s5m",
        "colab_type": "code",
        "colab": {}
      },
      "cell_type": "code",
      "source": [
        "df_bio_read.shape"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "metadata": {
        "id": "YDfurrGCputy",
        "colab_type": "text"
      },
      "cell_type": "markdown",
      "source": [
        "## **Verify the processing procedure (only if *Process \"raw\" data* was run)**"
      ]
    },
    {
      "metadata": {
        "id": "ZZD3V6zHnmvL",
        "colab_type": "code",
        "colab": {}
      },
      "cell_type": "code",
      "source": [
        "v = np.array(df_bio['log2_cleaned_ratio']).flatten() - df_bio_read['log2_cleaned_ratio']"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "metadata": {
        "id": "NpPSFdEh2JP3",
        "colab_type": "code",
        "colab": {}
      },
      "cell_type": "code",
      "source": [
        "min(v), max(v)"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "metadata": {
        "id": "pYgPlN60CMut",
        "colab_type": "code",
        "colab": {}
      },
      "cell_type": "code",
      "source": [
        "v = np.array(df_bio['log2_cleaned_ratio_th2d']).flatten() - df_bio_read['log2_cleaned_ratio_th2d']"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "metadata": {
        "id": "mF1KY5GSYvki",
        "colab_type": "code",
        "colab": {}
      },
      "cell_type": "code",
      "source": [
        "min(v), max(v)"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "metadata": {
        "id": "9zyQWSQnKutF",
        "colab_type": "code",
        "colab": {}
      },
      "cell_type": "code",
      "source": [
        "del df_bio"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "metadata": {
        "id": "P7qH-Zgbi7rG",
        "colab_type": "text"
      },
      "cell_type": "markdown",
      "source": [
        "## **Convert processed table back to design matrix format.**"
      ]
    },
    {
      "metadata": {
        "id": "X_ysPn7rNAQy",
        "colab_type": "code",
        "colab": {}
      },
      "cell_type": "code",
      "source": [
        "def design_from_biologist(df):\n",
        "  df.set_index(DESIGN_IDX+['GeneName'], inplace=True)\n",
        "  log2_columns = ['log2_' + x for x in LOG_SUFFIXES]\n",
        "  df[log2_columns] *= np.log(2.0)\n",
        "  cols = [x.replace('log2', 'log') for x in df.columns]\n",
        "  df.columns = cols\n",
        "  return df.unstack()"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "metadata": {
        "id": "HnUcyMkkKYj1",
        "colab_type": "code",
        "colab": {}
      },
      "cell_type": "code",
      "source": [
        "df_dm_read = design_from_biologist(df_bio_read)"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "metadata": {
        "id": "aEOgq4_SQEpS",
        "colab_type": "code",
        "colab": {}
      },
      "cell_type": "code",
      "source": [
        "df_dm = df_dm_read"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "metadata": {
        "id": "OPsRxkPPjMO0",
        "colab_type": "code",
        "colab": {}
      },
      "cell_type": "code",
      "source": [
        "del df_bio_read"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "metadata": {
        "id": "_URZORnmjO5q",
        "colab_type": "text"
      },
      "cell_type": "markdown",
      "source": [
        "## **Add expression level variables (not logarithmic).**"
      ]
    },
    {
      "metadata": {
        "id": "tYxHNKLJ7edf",
        "colab_type": "code",
        "colab": {}
      },
      "cell_type": "code",
      "source": [
        "def add_linear_column_blocks(df):\n",
        "  genes = df['r_g_median'].columns\n",
        "  for block in df.columns.levels[0]:\n",
        "    if block.startswith('log_'):\n",
        "      df = add_gene_block(df, np.exp(df[block]), block[4:], genes)\n",
        "  return df"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "metadata": {
        "id": "Ut1vEZ4M7hdr",
        "colab_type": "code",
        "colab": {}
      },
      "cell_type": "code",
      "source": [
        "df_dm = add_linear_column_blocks(df_dm)"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "metadata": {
        "id": "XmZ1M77orzgE",
        "colab_type": "text"
      },
      "cell_type": "markdown",
      "source": [
        "# **Linear modeling**"
      ]
    },
    {
      "metadata": {
        "id": "FDnXgMFUkACO",
        "colab_type": "text"
      },
      "cell_type": "markdown",
      "source": [
        "## **Construct the design for modeling.**"
      ]
    },
    {
      "metadata": {
        "id": "QFfKzn6KJ1D5",
        "colab_type": "code",
        "colab": {}
      },
      "cell_type": "code",
      "source": [
        "def construct_coefficient_labels(df, model_quadratic=True,\n",
        "                                 corrected_intercept=False):\n",
        "  \"\"\"Compute the column names for the design matrix.\"\"\"\n",
        "  # These are the names of the genes.\n",
        "  gene_labels = df['r_g_median'].columns\n",
        "\n",
        "  # These are the names of the coefficients.\n",
        "  # The quadratic terms need to be added.\n",
        "  coef_labels = []\n",
        "  coef_labels.extend(gene_labels)\n",
        "  if model_quadratic:\n",
        "    coef_labels.extend(['quad_' + x for x in gene_labels])\n",
        "  if corrected_intercept:\n",
        "    coef_labels.append('corrected_intercept')\n",
        "  return pd.Series(coef_labels)"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "metadata": {
        "id": "uNDI9luK0UPj",
        "colab_type": "code",
        "colab": {}
      },
      "cell_type": "code",
      "source": [
        "def construct_full_design(df,\n",
        "                          in_sample_f,\n",
        "                          out_of_sample_f,\n",
        "                          gene,\n",
        "                          delta_model=False,\n",
        "                          thresholded='',\n",
        "                          model_quadratic_terms=True,\n",
        "                          log_weights=False,\n",
        "                          inverse_design=False,\n",
        "                          corrected_intercept=False,\n",
        "                          x0_op=None,\n",
        "                          diff_op=None,\n",
        "                          inv_diff_op=None):\n",
        "  \"\"\"Construct the design matrix and dependent variable.\n",
        "\n",
        "  For a given gene, the self promoted experiments are filtered out.  Then,\n",
        "  the design matrix for the quadratic model is constructed.\n",
        "\n",
        "  Args:\n",
        "    df: data frame with time points as rows, genes as columns.\n",
        "    in_sample_f: Build model on this filter\n",
        "    out_of_sample_f: Compute out of sample deviance on this filter\n",
        "    gene: gene to build model for.\n",
        "    delta_model: whether to subtract assumed steady state solution.\n",
        "    thresholded: suffix tag for thresholded data.\n",
        "    model_quadratic_terms: whether to include quadratic terms\n",
        "    log_weights: whether to apply weights to make log model\n",
        "    inverse_design: (boolean) whether the design matrix is inverse differenced.\n",
        "    corrected_intercept: (boolean) Fit intercept with \"integrate\" and \"log\",\n",
        "    x0_op: matrix to produce vector of time zero points for each timecourse\n",
        "    diff_op: matrix for differencing operator\n",
        "    inv_diff_op: matrix for inverse differencing operator\n",
        "  Returns:\n",
        "    design matrix, dependent variable, time=0 values,\n",
        "    weights, row_index, col_index, in_sample_f, out_of_sample_f\n",
        "  \"\"\"\n",
        "  notlast_f = df.index.get_level_values('rseq') != 0\n",
        "  notzero_f = df.index.get_level_values('time') != 0\n",
        "  row_index = df.index\n",
        "\n",
        "  col_index = pd.Series(\n",
        "      construct_coefficient_labels(\n",
        "          df, model_quadratic=model_quadratic_terms,\n",
        "          corrected_intercept=corrected_intercept))\n",
        "  print('Promoter experiment filtered, shape = %s' % (df.shape,))\n",
        "\n",
        "  # Select column sets according to input booleans.\n",
        "  val = 'cleaned_ratio' + thresholded\n",
        "  yval = val\n",
        "  # Construct the design matrix d_m.\n",
        "  #\n",
        "  # df has multi-indexed columns, where the final piece is \"gene\",\n",
        "  # e.g. df['ratio'] has one column per gene.\n",
        "  #\n",
        "  # Matrices, with rows (strain X time) and columns for genes:\n",
        "  #\n",
        "  # u_m is the matrix of r/g values, divided by the t0 ratio\n",
        "  # uu is a matrix of the same shape with the quadratic terms.\n",
        "  #\n",
        "  # u is the column of r/g values for the gene in question.\n",
        "  u_m = np.array(df[val])\n",
        "  if model_quadratic_terms:\n",
        "    u = np.array(df[val][gene]).reshape(-1, 1)\n",
        "    uu = u * u_m\n",
        "    d_m = np.hstack([u_m, uu])\n",
        "  else:\n",
        "    d_m = u_m\n",
        "\n",
        "  # We now subtract 1.0 from the design matrix, giving it\n",
        "  # column blocks like (u-1)(uu-1).\n",
        "  # The intercept is now the deviation from steady state at t=0.\n",
        "  # For a delta model we assume steady state, with intercept=0.\n",
        "  d_m -= 1.0\n",
        "\n",
        "  # We need to add a column for a non-standard intercept.\n",
        "  # This will be transformed by at least one of the log_weights and\n",
        "  # inverse_design clauses immediately below.\n",
        "  if corrected_intercept:\n",
        "    d_m = np.append(d_m, np.ones((d_m.shape[0], 1)), axis=1)\n",
        "\n",
        "  # For log models, we now divide by the gene in question.\n",
        "  if log_weights:\n",
        "    u = np.array(df[val][gene]).reshape(-1, 1)\n",
        "    d_m /= u\n",
        "    yval = 'log_cleaned_ratio' + thresholded\n",
        "\n",
        "  # For inverse models, we now \"integrate\" the design matrix.\n",
        "  if inverse_design:\n",
        "    d_m = np.dot(inv_diff_op, d_m)\n",
        "\n",
        "  if inverse_design:\n",
        "    yfull = df[yval][gene]\n",
        "    # The inverse difference operator returns 0 for all timeseries 0 points.\n",
        "    # We thus subtract out the zero values, and filter the 0 rows.\n",
        "    yzero = yfull - np.dot(x0_op, np.array(yfull))\n",
        "    y = yzero[notzero_f]\n",
        "    d_mat = d_m[notzero_f]\n",
        "    row_index = df.index[notzero_f]\n",
        "    weights = None\n",
        "    in_sample_f = in_sample_f[notzero_f]\n",
        "    out_of_sample_f = out_of_sample_f[notzero_f]\n",
        "  else:\n",
        "    # The final derivative is linearly dependent on the previous ones,\n",
        "    # because there are only N-1 differences.  Thus, we filter on rows that are\n",
        "    # not the final point in their timecourse.\n",
        "    y = np.dot(diff_op, df[yval][gene])[notlast_f]\n",
        "    d_mat = d_m[notlast_f]\n",
        "    row_index = df.index[notlast_f]\n",
        "    weights = 1.0 / (diff_op ** 2).sum(axis=1)[notlast_f]\n",
        "    in_sample_f = in_sample_f[notlast_f]\n",
        "    out_of_sample_f = out_of_sample_f[notlast_f]\n",
        "\n",
        "  y_base = np.dot(x0_op, np.array(df['r_g_median'][gene]))\n",
        "\n",
        "  return (d_mat, y, y_base, weights, row_index, col_index,\n",
        "          in_sample_f, out_of_sample_f)"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "metadata": {
        "id": "06ZDyYfs44vq",
        "colab_type": "code",
        "colab": {}
      },
      "cell_type": "code",
      "source": [
        "# Picking a gene to model - choose one that is a transcption factor\n",
        "gene = 'FKH1'\n",
        "promoted = df_dm.index.get_level_values('TF') == gene\n",
        "\n",
        "# Only hold out the self-experiment\n",
        "in_sample_f = np.full(len(df_dm.index), True, dtype=bool)\n",
        "out_of_sample_f = np.full(len(df_dm.index), False, dtype=bool)\n",
        "\n",
        "in_sample_f = in_sample_f[~promoted]\n",
        "\n",
        "out_of_sample_original_size_f = pd.Series(\n",
        "      out_of_sample_f & ~promoted, index=df_dm.index)\n",
        "out_of_sample_f = out_of_sample_f[~promoted]\n",
        "\n",
        "df = df_dm[~promoted]\n",
        "\n",
        "# Build the differencing operators\n",
        "times = np.array(df.index.get_level_values('time'))\n",
        "x0_op, diff_op, inv_diff_op = construct_difference_operators(times)"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "metadata": {
        "id": "sMNTw7Kz5t3C",
        "colab_type": "code",
        "colab": {}
      },
      "cell_type": "code",
      "source": [
        "(design_m, y_dep, y_base, weights, row_index,\n",
        " col_index, in_sample_f, out_of_sample_f) = construct_full_design(\n",
        "       df,\n",
        "       in_sample_f,\n",
        "       out_of_sample_f,\n",
        "       gene,\n",
        "       delta_model=True,\n",
        "       thresholded='_zth2dfilt',\n",
        "       model_quadratic_terms=True,\n",
        "       log_weights=True,\n",
        "       inverse_design=False,\n",
        "       corrected_intercept=False,\n",
        "       x0_op=x0_op,\n",
        "       diff_op=diff_op,\n",
        "       inv_diff_op=inv_diff_op)"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "metadata": {
        "id": "7o44WnWD6wg6",
        "colab_type": "code",
        "colab": {}
      },
      "cell_type": "code",
      "source": [
        "design_m.shape"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "metadata": {
        "id": "2nLBZpOBH1vG",
        "colab_type": "code",
        "colab": {}
      },
      "cell_type": "code",
      "source": [
        "y_dep.shape"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "metadata": {
        "id": "tO6fCruNr9fy",
        "colab_type": "text"
      },
      "cell_type": "markdown",
      "source": [
        "## **Build a simple linear model (sklearn, not glmnet).**"
      ]
    },
    {
      "metadata": {
        "id": "FWrsAPDxama2",
        "colab_type": "code",
        "colab": {}
      },
      "cell_type": "code",
      "source": [
        "# We are attempting to match our weighting scheme with GLMnet, imperfectly.\n",
        "larsmodel = sklearn.linear_model.Lars(fit_intercept=False, normalize=False)"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "metadata": {
        "id": "UaMVDpp0biGe",
        "colab_type": "code",
        "colab": {}
      },
      "cell_type": "code",
      "source": [
        "larsmodel.fit(design_m, y_dep)"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "metadata": {
        "id": "mQ3zAIi63wQ3",
        "colab_type": "code",
        "colab": {}
      },
      "cell_type": "code",
      "source": [
        "lambdacol = 40"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "metadata": {
        "id": "5itfBzglcN-2",
        "colab_type": "code",
        "colab": {}
      },
      "cell_type": "code",
      "source": [
        "f = larsmodel.coef_path_[:, lambdacol] != 0"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "metadata": {
        "id": "DgdC0uaedT-e",
        "colab_type": "code",
        "colab": {}
      },
      "cell_type": "code",
      "source": [
        "coefs = (larsmodel.coef_path_)[f, lambdacol]"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "metadata": {
        "id": "wRSSfPSXd4t8",
        "colab_type": "code",
        "colab": {}
      },
      "cell_type": "code",
      "source": [
        "for gene, coef in zip(col_index[f], coefs):\n",
        "  print('%15s % .4e' % (gene, coef))"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "metadata": {
        "id": "mMUODykqsit6",
        "colab_type": "text"
      },
      "cell_type": "markdown",
      "source": [
        "# **Prediction model**"
      ]
    },
    {
      "metadata": {
        "id": "vnw0jeZ7jkSi",
        "colab_type": "text"
      },
      "cell_type": "markdown",
      "source": [
        "## **Read in the prediction model.**"
      ]
    },
    {
      "metadata": {
        "id": "yGkXFKsE_grA",
        "colab_type": "code",
        "colab": {}
      },
      "cell_type": "code",
      "source": [
        "# Read in the prediction model.\n",
        "\n",
        "# Before validation experiments.\n",
        "coeffs = pd.read_csv(os.path.join(datadir, 'insample_coefs_20170601.csv'))\n",
        "\n",
        "# After validation experiments.\n",
        "#coeffs = pd.read_csv(os.path.join(datadir, 'insample_coefs_20180826.csv'))"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "metadata": {
        "id": "TZ9TuFhss2_A",
        "colab_type": "text"
      },
      "cell_type": "markdown",
      "source": [
        "## **Construct predictions from model.**"
      ]
    },
    {
      "metadata": {
        "id": "L-CoXjDis4X2",
        "colab_type": "code",
        "colab": {}
      },
      "cell_type": "code",
      "source": [
        "design = df_dm['cleaned_ratio_zth2dfilt']"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "metadata": {
        "id": "ejVzTlrys-Ye",
        "colab_type": "code",
        "colab": {}
      },
      "cell_type": "code",
      "source": [
        "design.shape"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "metadata": {
        "id": "3U0KMD3XtkD2",
        "colab_type": "code",
        "colab": {}
      },
      "cell_type": "code",
      "source": [
        "x0op, diff, diff_inv = construct_difference_operators(np.array(design.index.get_level_values('time')))"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "metadata": {
        "id": "K8VU1sXSt1XD",
        "colab_type": "code",
        "colab": {}
      },
      "cell_type": "code",
      "source": [
        "# Predict for for each gene in turn.\n",
        "gene = 'YHB1'\n",
        "\n",
        "# Don't predict the TF in its own experiment!\n",
        "tf_filt = design.index.get_level_values('TF') != gene\n",
        "# Data for gene in question\n",
        "actualg = design[gene]\n",
        "actualg_diff = np.matmul(diff, np.log(actualg))\n",
        "# Consider coefficients affecting the gene.\n",
        "f = coeffs['effect'] == gene\n",
        "n_c = sum(f)\n",
        "# 2 extra columns for data, full prediction\n",
        "marginal_matrix = np.zeros((len(actualg), n_c + 2))\n",
        "# Fill in the data\n",
        "marginal_matrix[:, 0] = actualg_diff\n",
        "# Initialize the full prediction.\n",
        "prediction = 0. * actualg\n",
        "# Iterate over the coefficients that affect this gene.\n",
        "for i, x in enumerate(coeffs[f].iterrows()):\n",
        "  cause = x[1]['cause']\n",
        "  coef = x[1]['coef'] * tf_filt\n",
        "  # The selected model enforces intercept == 0, and is \"log\".\n",
        "  # Compute the marginal contribution to (d/dt) ln(gene)\n",
        "  if cause.startswith('quad'):\n",
        "    marginal = (design[cause[5:]] - 1.0 / design[gene]) * coef\n",
        "  else:\n",
        "    marginal = (design[cause] - 1.0) / design[gene] * coef\n",
        "  marginal_matrix[:, i + 2] = marginal\n",
        "  prediction += marginal\n",
        "# \"Integrate\" with the diff_inv operator to get ln(gene), then exponentiate\n",
        "data_pred = np.exp(np.matmul(diff_inv, prediction))\n",
        "# Fill in the predicted results\n",
        "marginal_matrix[:, 1] = prediction\n",
        "# Write out marginals\n",
        "columns = ['data', 'prediction'] + list(coeffs[f]['cause'])\n",
        "marginal_diff_df = pd.DataFrame(marginal_matrix, index=design.index, columns=columns)\n",
        "marginal_df = pd.DataFrame(np.exp(np.matmul(diff_inv, marginal_matrix)), index=design.index, columns=columns)"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "metadata": {
        "id": "kw84qgfQtY-o",
        "colab_type": "text"
      },
      "cell_type": "markdown",
      "source": [
        "## **Investigate individual timecourse.**"
      ]
    },
    {
      "metadata": {
        "id": "CFK-bIrIv2xM",
        "colab_type": "code",
        "colab": {}
      },
      "cell_type": "code",
      "source": [
        "aft1 = marginal_df.query('TF==\"AFT1\"')"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "metadata": {
        "id": "wgZmSOilwJ3M",
        "colab_type": "code",
        "colab": {}
      },
      "cell_type": "code",
      "source": [
        "hits = np.sum(np.abs(np.log(aft1)) > 0.2) > 0"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "metadata": {
        "id": "JHUoOajUt6Hq",
        "colab_type": "code",
        "colab": {}
      },
      "cell_type": "code",
      "source": [
        "aft1[aft1.columns[hits]]"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "metadata": {
        "id": "wP0TkyReuJ7Z",
        "colab_type": "code",
        "colab": {}
      },
      "cell_type": "code",
      "source": [
        ""
      ],
      "execution_count": 0,
      "outputs": []
    }
  ]
}
